#include <iostream>
#include <vector>
#include <cmath>
#include "Spline.H"
#include "SnippetsForAssignment.cpp"

int main(int argc, const char **argv)
{   
    using namespace std;
    using rVec = Vec<Real, 2>;
    using T_Poly = Polynomial<2, rVec>;
    std::cout << "------------------------Assignment B--------------------------" << endl;
    int B_PossibleN[5]={5,10,20,40,80};
    RealVect B_ErrorComplete(5),B_ErrorNotAKnot(5);
    for(int k=0;k<5;k++)
    {   
        int N = B_PossibleN[k];
        RealBox Box_complete,Box_notAknot;
        Box_complete.push_back({-1,funB(-1),funBder(-1)});
        Box_notAknot.push_back({-1,funB(-1)});
        for (int i = 1; i < N; i++)
        {
            Box_complete.push_back({-1 + 2.0 * i / N, funB(-1 + 2.0 * i / N)});
            Box_notAknot.push_back({-1 + 2.0 * i / N, funB(-1 + 2.0 * i / N)});
        }
        Box_complete.push_back({1,funB(1),funBder(1)});
        Box_notAknot.push_back({1,funB(1)});
        InterpConditions IC_c(Box_complete),IC_n(Box_notAknot);
        Spline<1,4,ppForm> spl1=interpolate<4>(IC_c,BCType::complete);
        Spline<1,4,ppForm> spl2=interpolate<4>(IC_n,BCType::notAknot);
        ///compute the error
        B_ErrorComplete[k]= 0,B_ErrorNotAKnot[k] = 0;
        for(int i=0; i<N ;i++)
        {
            B_ErrorComplete[k] = std::max(B_ErrorComplete[k],funB(-1+(2.0*i+1)/N)-spl1(-1+(2.0*i+1)/N)[0]);
            B_ErrorNotAKnot[k] = std::max(B_ErrorNotAKnot[k],funB(-1+(2.0*i+1)/N)-spl2(-1+(2.0*i+1)/N)[0]);
        };
        ///plot
        plotIn(spl1,"../../output/B_complete_" + std::to_string(k+1)+".m");
        plotIn(spl2,"../../output/B_notAknot_" + std::to_string(k+1)+".m");
    };
    ///generate the result
    std::cout << "For Complete situation, " << endl;
    displayErrorConvergence(B_ErrorComplete);
    std::cout << "For Not_a_knot situation, " << endl;
    displayErrorConvergence(B_ErrorNotAKnot);
    plotAgainstRunge_B("complete");
    plotAgainstRunge_B("notAknot");
    std::cout<<"Plots have been generated as '/Output/B_[XX]Plot.m'.\n";


    std::cout << "------------------------Assignment C--------------------------" << endl;
    RealBox Box1;
    InterpConditions IntpCondsBox1;
    Box1.push_back({-5,funC(-5)});
    for (int i = 1; i <= 10; i++)
        Box1.push_back({i-5.5, funC(i-5.5)});
    Box1.push_back({5,funC(5)});
    IntpCondsBox1.ResetIntpInfo(Box1);
    Spline<1,3,cardinalB> bspl1=interpolate<3>(IntpCondsBox1,middleSite,cardinalB);
    plotIn(bspl1,"../../output/C_459.m");
    RealBox Box2;
    InterpConditions IntpCondsBox2;
    Box2.push_back({-5,funC(-5),funCder(-5)});
    for (int i = 2; i <= 10; i++)
        Box2.push_back({-6.0+ i, funC(i-6)});
    Box2.push_back({5,funC(5),funCder(5)});
    IntpCondsBox2.ResetIntpInfo(Box2);
    Spline<1,4,cardinalB> bspl2=interpolate<4>(IntpCondsBox2,complete,cardinalB);
    plotIn(bspl2,"../../output/C_458.m");
    ///plot
    plotAgainstRunge_C("458");
    plotAgainstRunge_C("459");
    std::cout<<"Plots have been generated as '/Output/C_[XX].m'.\n";

    std::cout << "------------------------Assignment D--------------------------" << endl;
    Real C_Sites[7]={-3.5,-3,-0.5,0,0.5,3,3.5};
    ///compute the error
    for(int i=0;i<7;i++)
        std::cout<<"x = "<<C_Sites[i]
        <<": E_458 = "<< std::abs(bspl2(C_Sites[i])-funC(C_Sites[i]))
        <<", E_459 = "<< std::abs(bspl1(C_Sites[i])-funC(C_Sites[i]))<<'\n';


    std::cout << "------------------------Assignment E--------------------------" << endl;
    int E_PossibleN[3]={10,40,160};
    for(int k=0;k<3;k++)
    {
        vector<rVec> Box(E_PossibleN[k]+1);
        for (int i=0 ;i <= E_PossibleN[k];i++)
            Box[i]=heartCurve(2.0* i/E_PossibleN[k]);
        Spline<2, 4, ppForm> spl4 = fitCurve<4>(Box, BCType::periodic);
        plotIn(spl4,"../../output/E_heart_"+std::to_string(E_PossibleN[k])+".m");
        plotCharacteristics_E(Box,std::to_string(E_PossibleN[k]));
    }
    std::cout<<"Plots have been generated as '/Output/E_[XX].m'.\n";
    
    
    // const int N3 = 20;
    // RealBox Box1;
    // InterpConditions IntpCondsBox4;
    // for (int i = 0; i <= N3; i++)
    // {
    //     Box4.push_back({-1 + 2.0 * i / N3, funB(-1 + 2.0 * i / N3)});
    // }
    // IntpCondsBox4.ResetIntpInfo(Box4);
    // Spline<1, 2, ppForm> spl7 = interpolate<2>(IntpCondsBox4);
    // plotIn(spl7, "spl7.m", 500);



    // string bufB = "AssignmentB.m";
    // char file_nameB[20];
    // strcpy(file_nameB, bufB.c_str());
    // ofstream fout(file_nameB);
    // fout << "%This is a file about plotting AssignmentB." << endl;
    // fout << "figure(1);" << endl;
    // fout << "clear;x=-5:0.01:5;f=1./(1+x.*x);plot(x,f,'LineWidth',3);hold on" << endl;
    // for (int n = 0; n < 4; n++)
    // {
    //     fout << polyBox[n] << ";" << endl;
    //     fout << "plot(x,px,':','LineWidth',2);" << endl;
    // }
    // fout << "xlabel('x'),ylabel('y'),title('Runge Phenomenon');saveas(1,'AssignmentB.png');";
    // fout.close();
    // cout << "Generate AssignmentB.m!" << endl;
    // polyBox.clear();

    // cout << "-----------------------Assignment C-----------------------" << endl;
    // for (int n = 0; n < 4; n++)
    // {
    //     int N = 5 * (n + 1);
    //     RealBox tempBox;
    //     InterpConditions IntpCondsC;
    //     for (int i = 0; i < N; i++)
    //     {
    //         tempBox.push_back({ChebyshevZeros(N)[i], funC(ChebyshevZeros(N)[i])});
    //     }
    //     IntpCondsC.ResetIntpInfo(tempBox);
    //     Solver.ResetIntpConds(IntpCondsC);
    //     Solver.QuickInterp();
    //     polyBox.push_back(Solver.getInterPoly());
    // }
    // for (int n = 0; n < 4; n++)
    //     cout << polyBox[n] << endl;
    // string bufC = "AssignmentC.m";
    // char file_nameC[20];
    // strcpy(file_nameC, bufC.c_str());
    // ofstream fout1(file_nameC);
    // fout1 << "%This is a file about plotting AssignmentC." << endl;
    // fout1 << "figure(2);" << endl;
    // fout1 << "clear;x=-1:0.01:1;f=1./(1+25*x.*x);plot(x,f,'LineWidth',3);hold on" << endl;
    // for (int n = 0; n < 4; n++)
    // {
    //     fout1 << polyBox[n] << ";" << endl;
    //     fout1 << "plot(x,px,':','LineWidth',2);" << endl;
    // }
    // fout1 << "xlabel('x'),ylabel('y'),title('Chebyshev Interpolation');saveas(2,'AssignmentC.png');";
    // fout1.close();
    // cout << "Generate AssignmentC.m!" << endl;
    // polyBox.clear();

    // cout << "-----------------------Assignment D-----------------------" << endl;
    // InterpConditions IntpCondsD({{0,0,75},{3,225,77},{5,383,80},{8,623,74},{13,993,72}});
    // Solver.ResetIntpConds(IntpCondsD);
    // Solver.QuickInterp();
    // polyBox.push_back(Solver.getInterPoly());
    // cout<<polyBox[0]<<endl;
    // polyBox.push_back(polyBox[0].derivate());
    // cout<<"The derivation of "<<polyBox[1]<<endl;
    // cout<<"(a) When t = 10, "<<endl;
    // cout<<"Position: "<<polyBox[0].eval(10.0)<< " feet from Time 0."<<endl;
    // cout<<"Speed: "<<polyBox[1].eval(10.0)<< " feet per second."<<endl;
    // cout<<"(b) For t from 0 to 13, "<<endl;
    // bool is_overspeed = false;
    // for(int i=0; i<1300; i++)
    // {
    //     if(polyBox[1].eval(0.01*i)>81.01)
    //     {
    //         is_overspeed = true;
    //         cout << "the car exceeds the speed limit when t = "<<0.01*i<<" ."<<endl;
    //         break;
    //     };
    // }
    // if(is_overspeed == false)
    //     cout << "the car never exceeds the speed limit. Innocent! "<<endl;
    // polyBox.clear();

    // cout << "-----------------------Assignment E-----------------------" << endl;
    // InterpConditions IntpCondsE1({{0,6.67},{6,17.3},{10,42.7},{13,37.3},{17,30.1},{20,29.3},{28,28.7}});
    // InterpConditions IntpCondsE2({{0,6.67},{6,16.1},{10,18.9},{13,15.0},{17,10.6},{20,9.44},{28,8.89}});
    // Solver.ResetIntpConds(IntpCondsE1);
    // Solver.QuickInterp();
    // polyBox.push_back(Solver.getInterPoly());
    // Solver.ResetIntpConds(IntpCondsE2);
    // Solver.QuickInterp();
    // polyBox.push_back(Solver.getInterPoly());
    // for (int n = 0; n < 2; n++)
    //     cout << polyBox[n] << endl;
    // string bufE = "AssignmentE.m";
    // char file_nameE[20];
    // strcpy(file_nameE, bufE.c_str());
    // ofstream fout2(file_nameE);
    // fout2 << "%This is a file about plotting AssignmentE." << endl;
    // fout2 << "figure(3);" << endl;
    // fout2 << "clear;x=0:0.01:28;hold on;" << endl;
    // fout2 << "A"<< polyBox[0] << ";" << endl;
    // fout2 << "plot(x,Apx,'LineWidth',2);" << endl;
    // fout2 << "B"<< polyBox[1] << ";" << endl;
    // fout2 << "plot(x,Bpx,'LineWidth',2);" << endl;    
    // fout2 << "legend('Sp1','Sp2'),xlabel('x'),ylabel('y'),title('Average Weight Curves (0,28)');saveas(3,'AssignmentE1.png');"<<endl;
    // fout2 << "figure(4);" << endl;    
    // fout2 << "clear;x=28:0.01:43;hold on;" << endl;
    // fout2 << "A"<< polyBox[0] << ";" << endl;
    // fout2 << "plot(x,Apx,'LineWidth',2);" << endl;
    // fout2 << "B"<< polyBox[1] << ";" << endl;
    // fout2 << "plot(x,Bpx,'LineWidth',2);" << endl;    
    // fout2 << "legend('Sp1','Sp2'),xlabel('x'),ylabel('y'),title('Average Weight Curves (28,43)');saveas(4,'AssignmentE2.png');";
    // fout2.close();
    // cout << "Generate AssignmentE.m !" << endl;
    // polyBox.clear();
    // return 0;
}


